Diffusion limited aggregation as a Markovian process: 
site-sticking conditions 
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Cylindrical lattice diffusion limited aggregation (DLA) , with a narrow width N, is solved for site- 
sticking conditions using a Markovian matrix method (which was previously developed for the bond- 
sticking case). This matrix contains the probabilities that the front moves from one configuration 
to another at each growth step, calculated exactly by solving the Laplace equation and using the 
proper normalization. The method is applied for a series of approximations, which include only a 
finite number of rows near the front. The fractal dimensionality of the aggregate is extrapolated to 
a value near 1.68. 



I. INTRODUCTION 

Diffusion limited aggregation (DLA) |jj has been the 
subject of extensive study since it was first introduced. 
This model exhibits a growth process that produces 
highly ramified self similar patterns, which are believed 
to be fractals H . It seems that DLA captures the essen- 
tial mechanism in many natural growth processes, such 
as viscous fingering ||], dielectric breakdown j|]JJ, etc. 
In spite of the apparent simplicity of the model, an an- 
alytic solution is still unavailable. Particularly, the ex- 
act value of the fractal dimension is not known. Some 
of the analytic approaches employed so far include the 
fixed scale transformation (FST) [pi, real space renor- 
malizationgroup (RSRG) J7|,|||9| and conformal mapping 
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In DLA there is a seed cluster of particles fixed some- 
where. A particle is released at a distance from the clus- 
ter, and performs a random walk until it attempts to pen- 
etrate the fixed cluster, in which case it sticks. Then the 
next particle is released and so on. There are two com- 
mon types of sticking conditions. The sticking condition 
described above is called "bond-DLA" , because it occurs 
when the random walker attempts to cross a perimeter 
bond between an unoccupied site and the aggregate. In 
an earlier paper [ft4j, we solved the bond-DLA problem 
using a Markovian process. Here we apply similar meth- 
ods to the "site-DLA" case, where sticking occurs as soon 
as the random walker arrives at a site that is a nearest 
neighbor to the aggregate. Since it is believed that the 
large scale structure of DLA is not sensitive to the type 
of sticking conditions used [p>p|] , one expects both prob- 
lems to yield the same asymptotic fractal structures. 

DLA can be grown in various geometries. In this paper 
we deal with the cylindrical geometry in two dimensions 
(2D), where the particles are emitted from a distant hor- 
izontal line at the top, while the seed cluster is a parallel 
line at the bottom, with periodic boundary conditions 
on the sides. We only consider relatively narrow cylin- 



ders, with widths ranging from N = 2 to N = 10. Even 
though the analysis in this paper is solely 2D, the same 
techniques can be applied in higher dimensions. 

An exact solution of bond-DLA with N = 2 was pub- 
lished in 1998 PH . A generalization of the same approach 
was used in order to solve slightly wider cases with N be- 
tween 3 and 7 |0|. The solution presented in the latter 
case is not exact, but still, it presents a well controlled 
series of approximations in the sense that any desired 
numerical accuracy could be obtained, provided that a 
sufficiently high-order of approximation is used. The dif- 
ficulty with performing a high-order calculation is that 
its complexity grows exponentially. 

The main idea in these references and in this paper is 
to follow the dynamics of the growing front. The shape 
of the interface determines the unique solution to the 
Laplace equation that determines the growth probabili- 
ties. The structure of the aggregate behind the interface 
is irrelevant and so is the history that led to the current 
interface. Each growth process changes the interface. We 
can therefore describe DLA as a Markovian flow in the 
space of interface configurations. The Markov states are 
the possible shapes of the interface, which are indexed 
by an integer, usually denoted by i or j. Pi(t + 1), the 
probability that the interface is in state i at time t + 1, 
depends only on the state of the interface at time t. The 
conditional transition probabilities from state j to state 
i make up the evolution matrix Eij, which is time in- 
dependent. Thus, the dynamics of the Markov chain is 
described by the Master equations, 



Pi(t- 



l) = E^i P i(*)' 



t = l,2, 



(1.1) 



Each matrix element E^j corresponds to a particular 
growth process, and the sum of j runs over all the in- 
terface configurations (whose number may be infinite). 
In order to fully describe the dynamics, it is necessary 
to calculate the probabilities of all the possible growth 
processes, for each of the possible initial configurations 
of the interface. We calculate the growth probabilities 



by solving the discrete Laplace equation on a lattice for 
a function <f>, which corresponds to the average density 
of random walkers, 

V 2 $ = 0. (1.2) 

In the dielectric breakdown model (DBM) |||| <& has an 
electrostatic meaning, so it is also commonly referred to 
as the "potential" . 



Usually, the equation set (1.1) is infinite because the 
number of possible shapes the interface may assume is 
unlimited. This may pose a problem for two reasons. 
For one, it is difficult to include all of the possibilities 
systematically. The case of bond-DLA with N = 2 is 
a counter-example, where the complete set of possible 
configurations can be easily characterized using a single 
parameter. This is because the interface has the shape 
of a step whose height j can be any nonnegative inte- 
ger pq . For N > 2, however, it is very difficult to 
parametrise the shape of the interface, even with the use 
of more than one parameter, because complex overhangs 
may occur |14| . The second problem is that even if it was 
possible to account for a complete infinite set of configu- 
rations, it would still be awkward to analyze the Markov 
process, e.g., finding its fixed point. Instead of account- 
ing for all the configurations, we make an approximation 
by employing some consistent truncation scheme on the 
list of configurations. In the O'th-order approximation 
we include only the top O rows of a configuration and 
truncate the rest; The list of configurations is sorted ac- 
cording to the maximal height difference, Am, between 
the lowest and highest particles on the interface pf. In 
the O'th-order approximation, only a finite set of con- 
figurations with Am < O are taken into account. The 
configurations with Am > O are truncated so that only 
their top O rows are taken into account (below the O'th 
row all the sites are considered to be occupied). This 
truncation does not have a noticeable effect on the up- 
ward growth probability (the growth probability at the 
tip), because of the exponential decay of the potential 
inside deep fjords. Because of this exponential decay, 
the approximation converges very fast as a function of 
O. Unfortunately, the number of configurations diverges 
exponentially with O, so that the calculation can be car- 
ried out only for relatively low-order (depending on the 
width N and on the strength of the computer) . 

In the case of site-DLA, the situation is a bit simpler 
than in bond-DLA, because it is generally harder for the 
random walker to penetrate deep into a fjord. A particle 
will only be able to enter fjords that are three sites wide 
or more, unlike the case of bond-DLA, where a particle 
can go into a single column fjord. This makes the solu- 
tion of site-DLA with N = 2 and N = 3 much simpler, 
because they both have only a finite number of interface 
configurations. The narrowest cylinder that can have an 
arbitrarily deep fjord (a configuration with an arbitrarily 
large Am) arises for N = 4, and thus there is an infinite 
number of configurations. However, there can be no fluc- 
tuations in the width of the fjord, so in this sense this case 



resembles the N = 2 case in bond-DLA. For N > 4 the 
approximation method must be used, but generally, for 
the same N and O the number of configurations in site- 
DLA is much smaller than in bond-DLA, so it is possible 
to perform higher-order calculations for wider cylinders. 

Once an order of approximation O is chosen, there 
is only a finite number of configurations, N c (N,0), 
which depends both on N and on O. The Markov pro- 
cess is then closed and irreducible. Closed means that 
Si=i Ei t j = 1 for j = 1, . . . , N c , and irreducible means 
that there is a finite probability to go from any initial 
state j to any final state i during a finite number of time 
steps. A basic theorem in Markov theory states that 
a closed and irreducible process necessarily has a single 
fixed point ||17). This fixed point represents the steady- 
state probabilities for the various interface configurations 
in the asymptotic time limit. The theorem is also true 
for an infinite number of states, so one can conclude that 
the unapproximated process also converges to a steady 
state. As mentioned, this steady state is characterized 
for example by a time independent average density p. 

The fixed point equations are, 



p: = Y, E ^ p h * = i,2,... 



(1.3) 



This means that P* is the normalized eigenvector of the 
evolution matrix E with an eigenvalue of 1. Once the 
steady state P* is calculated, it is possible to evaluate 
the steady-state average upward growth probability, 



N a 



(Pu P >* = ]T p jWi)> 



3=1 



(1.4) 



where Pup(j) is the total upward growth probability for 
configuration j. The average steady-state density of the 
aggregate is then given by 



p(N) = 



1 



N( Puv 



(1.5) 



Here, the density is written explicitly as a function of N. 
By p(N) we denote the true value of the density, in the 
limit O — > oo. We denote the result of the O'th-order 
approximation by p c (N, O). 

As mentioned, the number of configurations grows ex- 
ponentially with O and N, so it becomes infeasible to 
make the calculation for high values of O and N. We per- 
form calculations for N < 10. The calculated densities 
and the number of configurations are presented in Sec. 
IJ. We find that in order to obtain a relative accuracy of 
about 1CP 4 , it is necessary to go up to O = N — 2, or 
O = N - 1. This is achieved for N < 8, but for N = 9, 10 
it is too heavy a task for our computer resources. In spite 
of this, we are able to successfully extrapolate O to infin- 
ity for N = 9 and N = 10. The calculated densities are 
compared to direct measurements from cylindrical site- 
DLA simulations, and are found to be the same up to 



the accuracy of the simulation, which is about 0.01% or 
better. 

The fractal dimension D is extracted from the assump- 
tion that p(N) oc N~( d - D \ where d = 2 is the Euclidean 
dimension. In general one should also expect some cor- 
rections to scaling, especially for low TV's, i.e., 

(1.6) 



p(N) = AN- {d - D ^ (1 + B/N e + ...) 



where A and B are some constants, 9 represents the 
leading correction exponent, and the dots stand for a 
series of higher powers of 1/N . This scaling hypothesis 
is validated by both the analytic enumeration computa- 
tion and by the simulations, which were conducted up to 
N = 128. The best fit of such a model to the enumera- 
tion data results in an estimate of the fractal dimension 
of cylindrical DLA in 2D, D = 1.68 ± 0.01. The same 
fit is also performed with the simulation data, yielding 
D = 1.671±0.001. This is to be compared with the value 
D sw 1.66 often found in the literature ||. 

The differences between bond-DLA and site-DLA are 
manifested in the boundary conditions for the Laplace 
equation, and in the way the growth probabilities are ex- 
tracted from the potential $. The boundary conditions 
at the top are, 



<9$(m, n) 
lim y —^—t = 1, n = 0, 



,iV-l, 



(1.7) 



where m is the vertical direction (the growth direction), 
and fi denotes the periodic lateral direction. This de- 
scribes a uniform flux of incoming particles. In the origi- 
nal DBM papers Q|| a uniform potential is used instead 
of a uniform gradient, but if the distant boundary is very 
far, then the differences between the solutions for the two 
cases is exponentially small [ p.6| . The determination of 
the boundary conditions on the aggregate should be done 
with care. In the case of bond-DLA the potential is set 
to on the aggregate itself, while in site-DLA the poten- 
tial should be set to on nearest neighbors sites, i.e., on 
sites where growth might occur. Also, the derivation of 
the growth probabilities from the potential in site-DLA 
is done a bit differently than in bond-DLA, as explained 
in the next section. 

The Laplace equation with these boundary conditions 
can be solved exactly pip"*!. The idea is to divide the 
plane (or space in higher dimensions) into two parts: the 
upper part is an empty semi-infinite rectangle that be- 
gins at the row of the highest site on the lower boundary 
and continues upward ad infinitum. The lower boundary 
is the set of sites on which the potential $ is set to 0, 
and it depends on the type of sticking-conditions used: 
In the case of bond-sticking conditions the boundary is 
the aggregate itself, whereas in the case of site-sticking 
conditions it is the set of sites that are nearest neighbors 
to the aggregate. The lower part contains the aggregate 
and extends from the highest row downwards. The row 
that contains the highest particle in the aggregate is usu- 
ally set as a reference row with m — 0. Thus, in bond- 
DLA, the upper part has m > 0, and the lower part has 



m < 0, and in site-DLA, the upper part has m > 1, and 
the lower part has m < 1 , as explained in more detail 
in Sec. 0. Note that in either case the dividing row is 
considered to belong to both parts. In the upper part, 
it is possible to express the potentials in row m + 1 as a 
linear combination of those in row m, 

N-l 

$(m + l,ft) = 1+ ^ $(m, n')g N {n - ri). (1.8) 

n'=0 

This is especially useful for the bottom row of the upper 
part, m = or m = 1 (depending on the type of stick- 
ing conditions). The boundary Green's function g^fo), 
appearing in Eq. (|l.8[), is given by 



(1.9) 



1 N ~ X 
9N(n) = — Y^ e K ' cos ( k i n )- 

1=0 



The finite set of allowed wave- vectors ki = ^-l for 
I = 0, . . . , N — 1, is imposed by the horizontal periodic- 
ity. The factor m is related to ki through the dispersion 
relation 



sinh (k/2) = ± sin (fc/2) . 
or more explicitly, 



(1.10) 



e -«(fc) = 2 _ cos( - fc ) _ y^ 2 _ cos ( fc ))2 _ L ( ln ) 

An interesting property of the Green's function is that 

N-l 



J2 g N (n) = 1. 



(1.12) 



This property was proved algebraically in Ref. |16[ and 
was used in Refs. fl6||L4|] to check the computations of 
the Green's func tion. It is also used in the sample cal- 
culation of Sec. II B in the current paper for the same 
purpose. 

Usually, there is no general derivation for the solution 
in the lower part. In spite of that, the number of sites in 
the lower part is finite and not too large, so it is possible 
to simply write the equations for each of the potentials. 
The solution of the resulting finite and linear set of equa- 
tions is then straightforward. 

The paper is organized as follows: In Sec. [Iljwe present 
in detail the differences in the computation of the growth 
probabilities between bond and site sticking conditions. 
This presentation also explains the connection with the 
Laplace equation more rigorously. After that we perform 
a few sample calculations, for N = 2 and N = 3, in order 
to demonstrate the method presented in the introduc- 
tion. We then report the results of the computations for 
N between 4 and 10 for various orders of approximations 
O. We point out that the results collapse onto a uni- 
versal function that enables the extrapolation O — ► oo 



for N — 8, 9 and 10. This extrapolation reduces the 
error appreciably. In Sec. HI we present the simula- 



tion we made in order to verify our theoretical predic- 
tions. This presentation also explains how the boundary 
Green's function gAr(n) is used in some way as a prob- 
ability function, in order to make the simulation more 
efficient. Wc summarize in Sec. |V|. 



II. ENUMERATION 

Our computation method is referred to as enumera- 
tion, because it involves a systematic processing of some 
complete lists of configurations. 



A. The differences between site and bond sticking 
conditions 



Before proceeding with the actual calculations, we 
point out the differences in the computation of the growth 
probabilities between bond and site sticking conditions, 
because these differences are the essence of this paper. 
In order to point out the differences, we first review 
the method for computing the growth probabilities with 
bond-sticking con ditio ns. The first step is to solve the 
Laplace equation (1.2) on a lattice, where the boundary 
conditions are $ = on the aggregate, and d&/dm = 1 
on the distant boundary. In the case of cylindrical geom- 
etry there are periodic boundary conditions on the sides, 
see Fig. H]. The sticking probability per bond is then 
given by 



$>, 



Pb 



Ef/<V 



(2.1) 



where the subscript b refers to a perimeter bond and 
$b refers to the potential difference across such a bond. 
Because the potential is null on the aggregate, this dif- 
ference is equal to the value of the potential in a nearest 
neighbor site. Finally, the growth probability per site 
is computed by multiplying the sticking probability per 
bond by the number of bonds associated with the site, 
N b , see Fig. § 

Why does this procedure give the exact growth proba- 
bilities? In order to answer this question we must return 
to the original definition of DLA, that involves a single 
random walker. The random walker is injected into a ran- 
dom site near the remote boundary and it diffuses until it 
attempts to penetrate the aggregate, in which case it gets 
stuck. By "penetrate" we mean that it is not sufficient for 
the random walker to be in a nearest neighbor site to the 
aggregate, but that it has to attempt to cross a perimeter 
bond in order for it to stick. A possible way of measuring 
the growth probabilities for a particular interface config- 
uration is to send many random walkers, one after the 
other, and remove them after they stick. One has to keep 



track of how many particles get stuck in each site. Even- 
tually, the growth probabilities per site are estimated by 
the fraction of particles that got stuck in each site. In- 
stead of releasing the random walkers one at a time, it 
is more efficient to release many of them simultaneously, 
and let them perform a random walk without interacting 
with each other. Moreover, instead of releasing a large 
amount of particles in a single batch and waiting until 
all of them stick, it is also possible to inject them at a 
constant rate near the boundary, i.e., in each time step in- 
ject a new particle into each site near the boundary with 
a uniform probability r. The advantage in this way of 
performing the measurement is that after an initial equi- 
libration time the system arrives at a steady state, which 
is characterized by a time-independent average number 
of random walkers in all of the sites, including sites that 
are not near any of the boundaries. In the steady state 
the average number of random walkers entering into the 
system in each time step at the upper boundary is equal 
to the average number of random walkers vanishing out 
of the lower boundary. 

Denote the average number of random walkers in each 
site in the steady-state by <I>(m,n). The crucial point is 
that $ is time independent. In order to calculate $ we 
note that it satisfies the discrete Laplace equation 

= V 2 $(TO,n) = -4$(m,n)+ (2.2) 

$(m + l,ri) + $(m — l,n) 4- $(m, n + 1) + &(m,n - 1), 

because every random walker is equally probable to go to 
any one of its nearest neighbor sites. Thus, on a general 
lattice (or graph) the Laplace equation states that the 
value of $ at each site is equal to the mean value of <& 
on its nearest neighbor sites. Special care should be give 
to sites near the boundaries. Near the upper boundary 
each site has only three nearest neighbors, and particles 
are added at a constant rate r, therefore, 

$(m, n) — j [<I>(m, n - 1) + $(m, n + 1) 

+$(m- l,n) + $(m, n)] +r . (2.3) 

Note that the last term on the left before r is §>(m,n), 
instead of <&(to + l,n), because the particles that ran- 
domly choose to go up are unable to do so because of the 
boundary, and therefore they remain in the same place. 
Now, let us define <&(to + 1, n) = $(m, n) + 4r as a fic- 
titio us d ensity above the boundary. Then we see that 
Eq. (2J3) turns into the standard Laplace equation. This 
shows that instead of using the injection rate parameter 
r, it is possible to use the regular Laplace equation with 
the Neumann type boundary conditions that require the 
specification of the electric field, which corresponds to 
the difference in the potential across the upper bound- 
ary. Since the value of r does not change the growth 
probabilities, we are free to choose any value for it. If, 
for example, we choose r = 1/4 then the boundary con- 
ditions at the top are 

— = $(m+l,n)-$(m,n) = 1, (2.4) 

am 



for n — 0, 1, . . . , N — 1. We choose the upper bound- 
ary to be very far away from the lower one, because this 
simplifies the analytic expressions involved in the solu- 
tion of the Laplace equation, while leaving the sticking 
probabilities practically unchanged. 

Near the bottom boundary the situation is a bit dif- 
ferent. Each random walker that attempts to go into the 
aggregate is taken out of the system. The steady state 
equation for these sites is therefore 



$(m,n) = -^$(m',0, 



(2.5) 



where the sum is taken over all the sites (m', n') that are 
nearest neighbors (nn) to (m, n). At the lower boundary 
the situation is similar, because once again, we obtain 
the regular Laplace equation, if we choose the boundary 
conditions $ = on the aggregate itself. 

The growth probability in each site is evaluated as the 
average number of random walkers that stick in that site 
per unit time, normalized by the total number of par- 
ticles sticking in a time step across the total length of 
the lower boundary. An average of 1/4 of the particles 
in a site choose to go in each direction. Particularly, a 
fraction of 1/4 of the particles vanish after choosing to 
go via bonds that connect to the aggregate. The aver- 
age total number of particles sticking in a site would be 
a sum over all of its interface bonds, ^6 ^/4 — A r h$/4. 
In the steady-state, the average total number of sticking 
random walkers is equal to the average total number of 
random walkers injected into the system. Near the upper 
boundary r = 1/4 random walkers are injected into each 
of the N sites. Therefore the normalization factor is N/4 
and the growth probability in each site is N^/N. In 
Rcf. pM we arrive at the same result using the discrete 
Gauss theorem. 

The situation in site-DLA is different in the boundary 
conditions used near the aggregate, and in the expression 
for the growth probabilities. Now, random walkers never 
arrive at sites that are nearest neighbors to the aggregate, 
because as soon as they do they get stuck and are imme- 
diately removed from the system. We therefore impose 
$ = not on the aggregate itself, but rather on its near- 



est neighbor sites, see Fig. 



In general, the boundary 



for the Laplace equation is obtained by coating the aggre- 
gate with a layer of circled sites, as shown in the figure. 
This differentiates between the boundary of the aggre- 
gate itself and the boundary for the Laplace equation, in 
the sense that it is possible that two different aggregates 
would have the same boundary for the Laplace equation, 
see Fig. [|. One can think that a random walker does not 
interact directly with the boundary of the aggregate, but 
rather, it interacts with the circled sites that make the 
boundary for the Laplace equation. This means that a 
random walker that obeys site-sticking conditions cannot 
be used as a probe in any way to determine which of the 
two aggregates in the figure are present. Consequently, 
any two different aggregates that have the same bound- 



ary for the Laplace equation must have the exact same 
set of growth probabilities and can be therefore consid- 
ered as equivalent. Thus, from now on when we refer 
to an interface configuration, we relate to the shape of 
the boundary for the Laplace equation. The probabil- 
ity to be in such a configuration is a sum over all the 
underlying aggregate configurations. Another effect of 
the transition to the Laplace boundary is the narrowing 
of fjords. The padding of the aggregate by circled sites 
causes all the fjords to be narrower by two sites. Thus, 
a random walker can only penetrate into aggregates that 
have branches that are at least three sites apart. 

As in the case of bond-DLA, the sticking probabilities 
are evaluated as fractions of the average number of ran- 
dom walkers that stick per unit time. Only now, we must 
sum over bonds that lead into the site, rather than out 
of it, as is the case in bond-DLA. The growth probability 
per site is therefore, 



Psite{ 



(m,n) = — ^$(m',0, 



(2.6) 



where p S it c (m,n) is the total sticking probability at the 
perimeter site (m,n), see Fig. |L Unlike the case of 
bond-sticking conditions, where a single potential deter- 
mines the sticking probability in a particular site, now 
the potentials in several different sites contribute. This 
is because in bond-DLA the random walker sticks be- 
fore it moves out of the site, whereas in site-DLA the 
random walker sticks after it moves into it. This differ- 
ence gives the upper most tip of the aggregate even a 
greater advantage relative to bond-DLA, because a sin- 
gle particle tip gathers contributions from three sides in 
site-DLA, whereas in bond-DLA the only contribution 
is from above. This comes in addition to the screening 
property of the Laplace equation (common to both types 
of sticking conditions), which causes the sticking prob- 
abilities at the lower parts of the interface to decrease 
exponentially. 



B. Exact solutions for N — 2, 3 

The best way of explaining the enumeration method 
is by showing some sample calculations in detail. We 
present here the two simplest cases, namely, N = 2 and 
N = 3. In these relatively simple cases there is only a 
finite number of configurations, so it is possible to get an 
exact solution with no need for approximations. 

For N = 2, the interface of the aggregate itself has an 
infinite number of possible configurations, because it has 
the shape of a step whose height j can be any nonnega- 
tive integer jl6| . However, in site-DLA there are only two 
distinct states: j = and j > 0. The case j = refers 
to a flat interface, i.e., the two columns have the same 
height, and a growth process will create a step with j = 1, 
with probability 1. For any step size j ' > there are only 
two sites where a random walker may stick: above the 



highest particle in the aggregate, or on its side, see Fig. 
y. There is no possibility for the random walker to pene- 
trate into a fjord in N = 2, because it is too narrow, and 
the particle would stick at its entrance. The two config- 
urations are indexed by i = 1 and i = 2 respectively and 
are shown in Fig. |6|. 

We now begin building the evolution matrix E, by find- 
ing the growth probabilities for each of the two configura- 
tions. As mentioned, configuration i = 1 turns into i = 2 
with probability 1, hence E\ t \ — and E 2 ,i = 1- ft is 
important to keep track of the total upward growth prob- 
ability for each configuration, p up (i), that corresponds to 
events in which a newly stuck particle is higher than all 
of the particles in the aggregate. In this case p up (l) = f • 

In order to solve for i = 2, we first hav e to compute 
the Green's function according to Eq. (1.9), which gives 



g 2 (Q) = 2 - V2 = 0.5858, 
g 2 (l) = a/2 - 1 = 0.4142. 



(2.7) 



We check our calculations by veri fying that #2(0) + 
32(1) = 1, as expected from Eq. ( 1.12| ). The poten- 
tial $ near the growth sites can be expressed i n te rms 
of the variable 2 = 4>(1,0) according to Eq. (1.8), as 
shown in Fig. pL We usually set the row containing the 
highest particle in the aggregate as the reference row, 
with m = 0. Thus, the row m = 1 always contains the 
highest circled site that belongs to the Laplace bound- 
ary. Each of the sites in row m = 1 contributes to the 
potentials in the sites in row m = 2. The weight of the 
contribution is equal to the value of the Green's function 
gnin), where n is the horizontal distance between the 
contributing site in row m = 1 and the evaluated site in 
row m = 2. In this simple case there is only one site 
with a nonzero potential, namely, $(1,0), which is yet 
unknown, and which we denote by 2. The site (1,1) on 
its side is nearest neighbor to the aggregate and there- 
fore we set $(1, 1) = 0. Thus, the potential of the sites 
in row m = 2 have only a contribution from x. More 
specifically, $(2, 0) = 1 +52(0)2 because it is right above 
x, and $(2, 1) = 1 + 52(1)2 because it is removed by one 
site. The potential $(2, 0) does not contribute to any 
growth process, but is important for solving for x. The 
variable x is found using its Laplace equation, 



This time there are three different bonds coming from 
two sites: there are two bonds coming from (1, 0) and an 
additional one coming from (2, 1). This upward growth 
results in the same configuration, so 



£2,2 = Pu P (2) = \ [2x + 1 + g 2 {l)x] = ^^ 



0.8536. 
(2.10) 



This concludes the calculation of all of the growth pro- 
cesses. The resulting evolution matrix is 



E 



0.1464 

1 0.8536 



(2-11) 



We verify that the matrix is properly normalized by not- 
ing that the sum of the terms in each of its columns is 
equal to 1, i.e., 



][>*,, = 1, j = 1,2. 



(2.12) 



i=l 



The general theorem mentioned in the introduction en- 
sures the existence of a single eigenvector with an eigen- 
value of 1, or in other words, a fixed point vector P* that 
satisfies, 



EP* 



(2.13) 



The f act that the process is closed is manifested in Eq. 
( |2. 12 ). The process is also irreducible because there is 
a finite probability to go from any initial state to any 
final state during a finite number of time steps. The fact 
that there is a single fixed point implies that starting 
from any initial state, the system will converge to the 
fixed point. This fixed point represents the asymptotic 
time probabilities for seeing either one of the two possible 
configurations. 

The eigenvalues are the roots of the characteristic poly- 
nomial, 



Ao = l, 

Ai = — — = -0.1464, 



(2-14) 



4.t = 1+32(0)2 

2- V2 



= 0.2929. 



(2.8) 



Growth in site (0, 0) results in the flat configuration i = 1. 
It can only occur via one bond from site (1,0), denoted 
by a bold double arrow (4) in Fig. ||. Hence, 



2 2- a/2 
£1,2 = - = —^- = 0.1464, 



(2.9) 



where the denominator comes from the normalization 
factor N = 2. Growth can also occur in site (1,1). 



and the normalized fixed-point vector P* is given by 

(2.15) 



P* = ^^ = 0.1277, 
P* = ±2±\ZI = 0.8723. 



The steady-state weights enable us to calculate the aver- 
age upward growth probability, 

(Pup)* - E ^PupW = 12 | 7 = °- 8723 > ( 2 - 16 ) 

2 — 1 

which is connected to the mean density, 



P(2) 



6- V2 



2 (p up ) 



0.5732. 



(2.17) 



It is also possible to calculate the rate of convergence 
to the steady state. In general, the rate of convergence is 
determined by the largest eigenvalue of E, other than 1. 
Suppose that at time t — the state of the system differs 
from the steady state P(0) ^ P*. The difference vector 



v(0) = P(0) - P* 



(2.18) 



belongs to the linear subspace of vectors V = {v| J2 v i — 
0}, because both P(0) and P* are normalized probability 
vectors and th us th e sum of their components is equal to 
1. Now, Eq. (|2.12| ) ensures that V is an eigen-subspace 
of E, and as such it must contain at least one eigen- 
vector. Since in this simple case the space of configura- 
tions is only two-dimensional, then V is one-dimensional, 
and v(0) is necessarily an eigenvector with the eigenvalue 
Ai = —0.1464. After t time steps the state of the system 
is 



P(t) = E*P(0) 



Alv(O). 



(2.19) 



Therefore, the deviation from the steady state decays ex- 
ponentially, 



P(i) 



Alv(O) = (-l)' C -*v(0), 



where 



log|Ai 



0.5584. 



(2.20) 



(2.21) 



This means that a single time step is practically suffi- 
cient to arrive at the steady state. All of these theoretical 
predictions agree with results obtained from numeric sim- 
ulations, up to the accuracy of the simulation, which is 
better than 10" 5 . This dynamics is actually exactly the 
same as the first-order approximation of the frustrated 
climber model in Ref. p4j , except that the analysis of the 
temporal convergence is a little bit more refined there. 

The solution of the case N = 3 is also relatively sim- 
ple, because again there is only a finite number of growth 
configurations. This is because the width of the widest 
possible fjord is two sites, which is still insufficient for a 
random walker to penetrate, i.e., a random walker sticks 
as soon as it enters into a fjord. The three possible con- 
figurations are indexed in Fig. R. These are the same as 
the three configurations of the first-order approximation 
for bond-DLA with N = 3 G3). As in the example of 
N = 2, we proceed to calculate the probabilities for ev- 
ery growth process in each of the configurations. Once 
again, we first calculate the Green's function, 



53(1) 



33(0) 

53(2) = 



. 6-y / 21 

1-93(0) 
2 



0.4725, 

x/21-3 _ 
6 



0.2638. 



(2.22) 



The first configuration, i = 1, grows with probabil- 
ity 1 into configuration i = 2. Thus, E 2 ,i = 1 and 



Ei t i = i?3.i = 0, and also p up (l) = 1. The potential 
diagram for i = 2 is shown in Fig. |3[ Because of symme- 
try it is possible to conclude that $(1,0) = $(1,2) = x. 
The Laplace equation for x is 



4ar = x + l + [g&(0)+S3(l)]x, 
9-V21 



10 



0.4417. 



(2.23) 



The sticking probability at (0, 0) is x/3, because there is 
a single connecting bond, and because the normalization 
factor is 1/3 for this case. The resulting configuration 
is i = 3, however, a sticking event at (0, 2) also leads to 



3, so that E 3 2 



2^#i = 0.2945. The other 



3* 15 

possibility is an upward growth at (2, 1), that results in 



Pur 



the initial configuration i = 2. Thus, E2.2 
1 - £3,2 = ^iip 1 = 0.7055, and E 1>2 = 0. 

The potential diagram for i — 3 is shown in Fig. 
The Laplace equation is 



4:e = 1 + <?3(0)x, 



=> x = 



6 



= 0.2835. 



(2.24) 



A sticking event in (0, 1) leads to i = 1, therefore 
Ei.j, = x/3. The other possible sticking events at (1,0) 
or (1, 2) involve upward growths, that result in i = 2, i.e., 

£2,3 = Pup (3) = 1 - x/3 = 2±^H = 0.9055. This com- 
pletes the calculation of all of the element of the evolution 
matrix: 



E 



0.0945 

1 0.7055 0.9055 
0.2945 



(2.25) 



The normalized fixed point of the matrix is 

P* = [ 0.0210, 0.7562, 0.2227]. (2.26) 

This enables the computation of the average upward 
growth probability, and of the average density: 



(Pu P )* = E J=i ^*Pup(j) = 0.756245, 



p - W^T 



= 0.440774. 



(2.27) 



The second largest eigenvalue determines the character- 
istic time constant of the exponential convergence to the 
steady state, 



1 



log|Ai 



= 0.56. 



(2.28) 



C. Approximations for N > 3 



The two examples of the previous sections, for N — 2 
and N = 3, are special because there is only a finite 
number of possible configurations; a random walker can- 
not enter a fjord whose width is less than three sites when 
using site sticking conditions. The case N = 4 is the nar- 
rowest cylinder that can have a fjord that is three sites 
wide. Since this fjord can be arbitrarily deep, there is 
an infinite number of configurations. In spite of that, 
every configuration that has a fjord, which is more than 
one site deep, is uniquely determined by its depth, i.e., 
there is only one configuration with Am = 2, a single 
configuration with Am = 3, and in general: a single con- 
figuration with a specific Am, if Am > 2. The unique 
configuration with Am = 2 is shown in Fig. HG, along 
with the single configuration with a specific Am that is 
larger than two. Other than that, there are four possible 
configurations with Am = 1, which are shown in Fig. |llj, 
and finally, the trivial flat configuration, with Am = 0. 

This case resembles bond-DLA with N = 2 [[16| , in the 
sense that in both cases there is an infinite number of 
configurations, but this infinity can be represented using 
a single parameter. In Ref. Eq] this parameter is called 
"the step size" and is denoted by j, but actually it is 
the same as Am. The case of site-DLA with N = 4 is a 
bit different, because there are four configurations with 
Am = 1 instead of one. There is also a resemblance 
between the solution of the Laplace equation for the 
two cases, because in both cases the Laplace equation is 
solved on a single column with zero boundary conditions 
on the sides. Thus, in both cases there is an exponential 
decay of the potential inside the fjord, which is governed 
by the multiplicative factor e~ Kf — 2 — \/3- This enables 
us to treat the current case in an analogous way to the 
previous one. This could have given us analytic expres- 
sions for the Markovian matrix Etj, i,j — 1,2,... ,oo, 
for the steady-state vector P* , i — 1, 2, . . . , oo, and for 
the distribution of gaps inside the aggregate. However, 
we omit the presentation of this calculation because it is 
not of main interest of this work, and so we treat the case 
of iV = 4 in the same way as N > 4. 

For N > 4 the boundary may be complex, and it can- 
not be easily characterized because the width of a fjord 
can fluctuate and overhangs may appear. We therefore 
use the approximation scheme described in the introduc- 
tion, which was also used for bond-DLA p4[ . The cal- 
culation procedure involves going over all the possible 
configurations to some order, and calculating their set of 
growth probabilities. It is feasible to perform this task 
manually when the number of configurations is relatively 
small, but as N and O increase, the number of configura- 
tions grows exponentially and it becomes impractical to 
do so. We use the same computer program that was used 
for bond-sticking conditions, after making the necessary 
adjustments due to the site-sticking conditions. Manual 
calculations may still be important as test cases to check 



the operation of the program. 

The program goes over all of the possible configura- 
tions systematically. It starts with the trivial flat config- 
uration (Am = 0), which is indexed by j = 1. This con- 
figuration has only one possible growth process, which 
occurs with probability 1, that turns the interface into 
configuration j = 2, which has a single bump. The pro- 
gram then continues to j = 2 and analyzes its growth 
probabilities. Every growth process changes the shape of 
the boundary. Each time a particle sticks in a certain site, 
the program has to identify the newly formed configura- 
tion. In order to do so, it marks all of the nearest neigh- 
bors of the newly attached particle, because new particles 
may stick there. The new configuration is searched for 
in the existing list of configurations, which were already 
analyzed by the program. If it does not exist then it 
is added at the end of the list. In either case the pro- 
gram identifies the index of the resultant configuration i. 
Now, if the index of the original configuration is j, then 
the growth probability is stored in the matrix element 

E i,j- 

A configuration is characterized using the set of sites 

that are connected to infinity because these are the sites 
that are accessible to the random walker. Of course, 
any site that is higher than the highest site on a certain 
boundary is connected to infinity. Hence, it is sufficient 
to specify only the set of sites that are not higher than 
the highest site (the region m < 1). A single growth 
process may cause a whole region of sites to disconnect 
from infinity, for instance by sealing off an entrance to 
a fjord. This means that it is not sufficient to mark the 
nearest neighbors of a newly attached particle, but that 
it is necessary to recheck the complete set of sites that are 
connected to infinity. We perform this by an algorithm 
that marks this set recursively. 

Special care has to be taken for upward growth pro- 
cesses, because they may cause Am to exceed O. In case 
this happens, the bottom row of the configuration is trun- 
cated. Finally, symmetry has to be taken into account. 
Rotations around the axis of the cylinder and reflections 
about any vertical axis do not change the growth prob- 
abilities or the steady-state weights, so the set of all of 
the symmetric configuration are represented by a single 
canonical choice. More specifically, a configuration is rep- 
resented by a binary word that consists oi N x O digits 
that correspond to the sites; The empty sites that are on 
the exterior are given a value of 1, and the rest of the 
sites are assigned with zeros. We choose the canonical 
form as the word that has the maximal numerical value. 

After the complete list of configurations is processed, 
the calculation of the evolution matrix is completed, and 
it is closed, i.e., £\ Eij — 1 for every j. Then the steady- 
state vector P* is calculated iteratively by applying the 
evolution matrix many times on some initial state vec- 
tor. This method is much faster than any of the standard 
techniques for solving a set of linear equations, especially 
when the number of variables is very large. The next 
step is to calculate the average upward growth proba- 



bility, according t o E g. ( |l.4| ), and the average density, 
according to Eq. (|l.5| ). Our computer resources enabled 
us to conduct the enumeration only up to a finite order 
O max that depends on N. As explained above, for N = 2 
and N = 3 there exists a finite number of configurations, 
and higher order approximations are irrelevant. One may 
be surprised that we are able to reach O = 8 for N = 6, 
but do not reach such a high order for N — 4 and N = 5, 
because for sure, there are less configurations in the same 
order of approximation for lower N's. The reason is that 
very good convergence is achieved already for O = N, 
so we had little to gain by going to much higher orders, 
and we stop at O = N + 2 for TV = 4 and N = 5. The 
calculated densities p c (^, O) are presented in Table |, to- 
gether with the number of configurations N c . The Table 
also presents the extrapolation and simulation results. 



D. The extrapolation of the order of approximation 
to infinity, O — > oo 

Very good accuracy (about 10 -4 ) is also obtained for 
N = 8, even-though O max = N — 2 = 6. However, for 
N > 9 the results are not very accurate, because the 
maximal available order is only O = 5 for N = 9,10 
and O = 4 for N = 11, 12. In spite of that, we are 
able to arrive at a more precise estimation for N = 9,10 
by extrapolating O — ► oo. The extrapolation does not 
improve the accuracy of the cases N = 11,12 to a 
satisfactory level. Our aim is to deduce the value of 
p(N) — limo^oo p c (N, O) from the limited range of avail- 
able values for O. We start by noting that our data 
practically reached asymptotia for N — 4, 5, 6. We de- 
tect that the differences, p c (0,N) — p c {0 + 1,N), de- 
cay exponentially and thus conclude that the function 
/ = In [p c (N, O)fp(N) — 1], is very close to being linear. 
Substituting the parameterization / = /3 — aO/N we are 
able to extract the three unknowns, a, /?, and p(N) using 
at least three data points. For N = 6 and O = 4,5, and 
6, we find that (3 = 0.03 and a = 12.31. The value of 
p(6) turns out to be very close to the highest available 
approximation p c (6, 8). 

Scaling theory would imply that, for large N and O, f 
should become a universal function, which depends only 
on the scaled ratio x = O/N (without an additional de- 
pendence on N). Following this expectation, we thus 
conjecture the general relation 



p c (N,0)=p(N) 



1 



of (O/N) 



(2.29) 



with f(x) ~ -12.3x, for N, O > 1. 

To test this conjecture, we estimated p(N), for N > 4, 



p(N) 



Pc{N,O max ) 

1 + e /(O max /JV) ' 



(2.30) 



We have then used this estimate to calculate 
p c {N, 0)/p{N) for O < Omax- The resulting values are 



shown in Fig. |2|, together with the line f(x) = —12.3a;. 
Clearly, all the values for O/N > 0.4 are consistent with 
our conjectured form for f(x). 



The values of p(N), as deduced using Eq. ( p.30| ), are 
listed in Table Q. Clearly, they all agree with the values 
from the simulations, except for small deviations that 
appear for N — 9 and 10. In the cases N = 11, 12 the 
deviations are relatively large, because O max is too small, 
and hence the extrapoltion results are not specified. 



E. An enumeration based estimate of the fractal 
dimension D 



In the previous section, we obtain very accurate esti- 
mates of the asymptotic (O — > oo) average steady-state 
densities p(N). In this section, we extrapolate the latter 
densities in the limit N — > oo, in order to find the fractal 
dimension D. Consider a N d segment in the steady state 
regime of growth. Assuming that the structure is a self 
similar fractal, which has no characteristic length scale 
other than N, we expect that the average mass of the seg- 
ment would be proportional to N D , and that the density 
would be proportional to N D ~ d . In principle how eve r, 
one expects some corrections to scaling as in Eq. (jTjj). 
Taking only the first correction term of that equation into 
account we get an approximation that depends on four 
parameters: D, A, B and 8: 



Pa {N) - AN 



D-d 



B 

w 



(2.31) 



where the subscript a denotes that this is an approxi- 
mation. Usin g the four data points with 7 < N < 10, 
a fit to Eq. ( |2.3lD yields D = 1.64, log(A) = -0.63, 
B = 1.31 and 9 = 1.48. The calculation of the pa- 
rameters can also be based on more than the minimal 
four points, using a least mean square error method. We 
choose to minimize the logarithmic (or relative) errors 
Ap/p rather than the errors in the densities Ap, because 
we find them to be more uniformly distributed. The re- 
sults of the fit using the six data points with 5 < N < 10 
yields D = 1.74±0.06, log(A) = -1.0±0.3, B = 1.5±0.6, 
and 9 = 0.80±0.13. The error estimates are evaluated us- 
ing a confidence level of 0.95. Since the fit yields a value 
for 9 that is close to 1, we also try a three parameter fit, 
fixing 9 = 1. Using the three rightmost data points, for 
N = 8,9 and 10, gives D = 1.68, log(A) = -0.799, and 
B = 1.16. Using more points with 5 < N < 10, we get 



D = 1.68 ±0.01 

log(A) = -0.784 ±0.016 

B = 1.12 ±0.05. 



(2.32) 



Finally, an alternative four parameter form, including 
only "analytic" corrections, is 



Pa (N) = AN 



D-d 



B 

N 



C_ 

N 2 



(2.33) 



This time the results for 7 < N < 10 are D = 1.65, 
log(^l) = -0.68, B = 0.55 and C = 1.10, and the 
least mean square calculation for 5 < N < 10 yields 
D = 1.70 ± 0.02, log(^) = -0.87 ± 0.08, B = 1.5 ± 0.4, 
and C = -0.5 ± 0.4. 

We thus conclude that the fractal dimension of cylin- 
drical DLA is D ~ 1.68 ± 0.01, close to the results of 
earlier numerical work |6]| . 



III. SIMULATION 

As mentioned, our analytical enumeration results are 
confirmed by simulations. In this section we describe how 
our simulations were conducted, with special attention to 
the boundary Green's function g^(n), which is given a 
new probabilistic meaning. We also discuss the accuracy 
of the results, and finally, we try to fit the results to some 
approximations as in the end of the previous section and 
obtain some more estimates of the fractal dimension. 

Our simulation is performed on a lattice, which is rep- 
resented by a 2D array variable. Each of the variables 
in the array can assume one of two possible values, 1 or 
0, that determine whether the relevant site is occupied 
by an aggregate particle or not, respectively. The size 
of the array is (14A^) x N, i.e., its width is N and it is 
composed of 14 blocks of N x N sites stacked one on top 
of the other. The number 14 is quite arbitrary and could 
be chosen differently. In principle, the lattice should be 
tall enough to allow the aggregate to arrive at a steady 
state, and also to allow a margin at the top, because the 
average density of the aggregate is lower near the growing 
front. Each time a new cluster is initialized, the lattice 
array is cleared so that all of its variables are set to 0, 
except for the bottom row, which is set to 1. This means 
that the initial shape of the aggregate is a horizontal line 
at the bottom of the lattice. A random walker is char- 
acterized by the coordinates (to, n) of its position. In 
each simulation step a direction is chosen randomly and 
the particle is advanced in that direction. If the parti- 
cle happens to go into a site that is nearest neighbor to 
the aggregate then it sticks, i.e., the value of the relevant 
lattice variable is updated from to 1. Then the next 
random walker is released, and so on. 



A. The role of the Green's function 

In principle, each new random walker should be re- 
leased far above the aggregate, near the upper distant 
boundary. In practice, nothing can happen to the ran- 
dom walker (it cannot stick) until it crosses the bold line 
in Fig. O This line is drawn between the highest row 
where a random walker can stick (to = 1) and the row 
above it (m = 2), and thus it differentiates between the 



active zone below the line, with to < 1, and the inactive 
zone above it, with to, > 1. The projection of the path 
of the random walker on the vertical axis (its to coordi- 
nate) is also a random walk, only in one dimension (ID). 
Usually in ID there is a probability of 1/2 to go up and 
the same probability to go down, but in our case, there 
is a probability of 1/4 to go in cither direction, and a 
probability of 1/2 to stay at the same row. Nevertheless, 
this motion is still equivalent to a random walk, however 
the effective time step is longer. A quality of ID ran- 
dom walks is that there is a probability of 1 to arrive at 
any site (no matter how far) within a finite time. There- 
fore, there is a probability 1 that eventually the random 
walker would cross the line from the inactive zone into 
the active zone. The random walker is equally proba- 
ble to cross this line at any of the N sites, so instead of 
waiting for a long time, it is more efficient to start the 
simulation by inserting the random walker in a random 
site just below the line in the active zone |l|] . 

But what happens if the path of the random walker 
happens to cross the line into the inactive zone? Once 
more we apply the same reasoning and claim that ulti- 
mately the random walker would re-cross the line down- 
wards at some point with probability 1 . Unlike the initial 
insertion, this time the distribution of the reentry point 
is not uniform. It is quite easy to see for example, that 
there is a greater chance for the particle to reenter at the 
exact same site from which it exited than for it to reenter 
at a site that is far away. Let us denote by ^(to,?^ n') 
the probability that if the particle is at some initial site 
(m,n) in the inactive zone (to > 1), it will cross the line 
for the first time at (to' = l,n'). In the next time step 
the random walker moves to one of its nearest neighbors 
with equal probability. Therefore ^(to, n; n') must be 
equal to the average of ^f on all the nearest neighbors. 
This implies that W satisfies the Laplace equation (in the 
coordinates to and n), 



V 2 *(m, n;ri) = 0. 



(3.1) 



The boundary conditions for V? at the lower boundary 
are 



^(m — l,n; n') 



1 , n = n' 
, otherwise 



(3.2) 



This is true because if the random walker is already in 
the row to' = 1 then it already passed the line between 
to = 1 and to = 2 and so it stops before it starts. The 
boundary conditions at the top are ^ = const., or equiv- 
alently 



lim — = 0. 

twoo dm 



(3.3) 



These are the exact same conditions satisfied by the 
Green's function |lq] , and so the theorem about the 
uniqueness of the solution of the Laplace equation with 
boundary condition assures that ^ is equal to the Green 



10 



function, and especially at the first row above the line 
m = 1, 



*(l,n;n') = g N (n - n'). 



(3.4) 



This means that each time the random walker attempts 
to cross the line to the inactive zone, it can be returned to 
the active zone immediately. The distance of the reentry 
point from the exit point should be chosen randomly from 
the distribution defined by the Green's function gN(n). 
This policy saves a lot of simulation time in comparison 
with the alternative option of letting the random walker 
wander freely until it finally sticks, or until it passes some 
arbitrary critical distance from the aggregate. We note 
i n pa ssing that the discussion in this section proves Eq. 
(1.12) in an alternative, probabilistic approach, simply 
due to the fact that </jv(n) is a probability function. 



B. Analyzing the statistics 

A single cluster is completed as soon as the first parti- 
cle sticks in the top row of the lattice. Then, the number 
of particles in each row is counted and stored in a (ID) 
array variable that represents the average density profile 
as a function of height. Then the lattice array is cleared 
and a new aggregate is started. In contrast, the density 
profile array is not cleared, and it accumulates data for 
each new cluster so that after many iterations it con- 
verges to the average density, when normalized by the 
number of iterations iVj. 

An example of a density profile is shown in Fig. U4, 
where N = 10 and JVj sw 2 x 10 7 . Three distinct regions 
are visible in the graph; On the left part there is a fast 
decay from an initial density of 1 to a plateau. These 
graphs always start from a density of 1, because the ini- 
tial conditions for growth are that the bottom row of the 
lattice is completely occupied. The decay to the plateau 
shows the convergence to the steady state stage of the 
growth. It seems that the steady state settles roughly at 
a height that is equal to the width, i.e., about 10. The 
middle section of the graph seems to be a flat plateau of 
constant density. In fact, there are small statistical fluc- 
tuations due to the randomness of the simulations, which 
are invisible because they are on the order of 10~ 5 . Fi- 
nally, near the right end of the graph there is a decay to 
0. This is because the density near the growing front is 
smaller than in the frozen part, because more particles 
are still expected to stick there and finally raise the den- 
sity to the steady state value. Naturally, the density is 
above the highest particle in the aggregate, which is 
always at the top row of the lattice, because the simula- 
tion always stops at that point. It seems that the width 
of the interface layer is also close to N. 

Only the middle section of the aggregate is taken into 
account in measuring the average steady state density 
p(N). As mentioned, our impression is that a margin 
of N sites from each side is enough, but we work with 



margins of 2N. We thus evaluate p(N) as the average 
of the density profile over the plateau area. A possible 
way of estimating the accuracy of this average is by tak- 
ing the standard deviation and normalizing by the square 
root of the number of rows that participate in the aver- 
aging. This is however, somewhat optimistic and would 
produce very low error estimates, because this calculation 
assumes statistical independence between adjacent rows, 
where in fact, there are significant correlations. The right 
factor to normalize by is therefore the square root of the 
effective number of independent rows. Since the aggre- 
gate is fractal, the only available length-scale is N, and 
therefore the correlation length £ should be proportional 
to N. We therefore conservatively guess that two rows 
that are N rows apart are independent, and hence esti- 
mate the accuracy by the standard deviation divided by 
\/T0, because there are 10 blocks of N x N is the steady 
state region. This error estimate is based on the expected 
dependence on the number of blocks, but there could be 
some numerical factor missing. 

An alternative way of measuring p(N) is by measuring 
(Pup)* directly and using Eq. (|L~5|). After the aggre- 
gate reaches a height of 27V we assume that it is in the 
steady state and we start gathering statistics. In particu- 
lar, we count the number of upward growth events, when 
the random walker sticks above all the particles in the 
aggregate. Our results show very good correspondence 
between the two different ways; The typical relative dif- 
ference is on the order of 10~ 6 . 

The simulations were carried out for the following val- 
ues of TV: 2,3,..., 12, 16, 24, 32, 48, 64, 96 and 128. We 
did not go beyond that because our computer resources 
did not suffice to iterate a large enough number of clus- 
ters to obtain a relative accuracy of about 10 -4 or better, 
as obtained for the other cases. 

We now proceed to fit the results for the 10 available 
data poin ts with 128 > N > 10 in a similar way to 
Sec. HE. The difference is that now we use the error 



estimates at to give weights to the different data points, 
because not all the accuracies are the same. This way 
the fit will allow greater residuals for data points with 
larger error estimates. Our first attem pt is to fit the 
four parameter approximation of Eq. (2.31). The re- 
sults are D = 1.673 ± 0.002, log (A) = -0.770 ± 0.0013, 
B = 1.03 ±0.06, and 6 = 0.96 ±0.06. The maximal rela- 
tive residual is 1.2 x 10~ 4 . Once more we set 6 — 1 and 
perform the fit for the remaining three parameters. The 
results are 



D = 1.671 ±0.001, 
log(A) = -0.762 ±0.003, 
B = 1.071 ±0.015. 



(3.5) 



The resulting error estimates seem a bit too optimistic, 
perhaps also becuase of the presence of some systematic 
errors that are not taken into account. The simulation 
results are shown in Fig. [l5| on log-log scales as plus 
signs, along with the latter three-parameter fit, shown 
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as a dashed line. The figure also shows the enumeration 
results as circles. Since the differences are hardly no- 
ticeable, we display the relative (logarithmic) residuals 
Vi = (p(Ni) - pa(Ni)) / p(Ni) separately in Fig. |l| on 
semi-log scales, in comparison with the relative error es- 
timates ±<Jj. The maximal relative residual is 1.3 x 1CP 4 . 
This is consistent with the order of magnitude of the es- 
timated a priori errors. 

A factor that indicates the compatibility between the 
a priori error estimates <x; and the a posteriori residuals 



Vi is, 



2 _ 



-E 



(3.6) 



where Na, the number of degrees of freedom, is equal 
to the number of data-points minus the number of un- 
known parameters. The value of \ 2 should be close to 
1. In the latter fit we get \ 2 = 0-9: whereas \ 2 = 0-3 m 
the former. The results of the fit imply that the three 
parameter approximation is sound. 

For the sake of comparison we also try to fit to the 
other test approximations that were introduced in the 
previous section. The best fit to the four parameter 
approximation in Eq. (|2.33|) is D = 1.6721 ± 0.0012, 
log(A) = -0.766 ± 0.006, B = 1.12 ± 0.07 and C = 
—0.28 ± 0.4. In this approximation, the residuals are 
not lowered drastically; the maximal relative residual is 
1.1 x 10~ 4 and x 2 — 0.2. The error estimate of fourth 
parameter C is much greater than the error estimates of 
the other parameters. The contribution of the term with 
the C parameter is on the same order of magnitude as 
the residuals, at least for the data point with large TV's. 
This implies that there may be significant contributions 
from the noise (the errors) in these data points to the 
parameter, and therefore its inclusion is redundant. 



IV. SUMMARY 

In this paper, we continue our endeavour to solve cylin- 
drical DLA analytically, i.e, to calculate the steady state 
average density p, as a function of the cylinder width TV, 
and to find the fractal dimension D. Unlike our previous 
work, which deals with bond-sticking conditions ]14J , this 
work solves for site-sticking conditions. The immediate 
problem in following our Markovian method is that, ex- 
cept for TV = 2, 3, there is usually an infinite number of 
configurations. The case TV = 4 has an infinite number 
of configurations, but is still relatively simple. The large 
variety of possible complex interface shapes for TV > 5 
prevents the inclusion of all the configurations and com- 
pels the use of an approximation scheme, in which only 
a finite number of rows O of the growing front near the 
tip are included. This approximation works because of 
the exponential decay of the Laplace potential <& inside 
deep fjords. The approximation leaves a finite number of 



configurations to work with, and thus the computational 
procedure can be completed. 

We find that this is a well controlled approximation, 
in the sense that any desired numerical accuracy can be 
achieved provided that a high enough order of approxi- 
mation O is used. The results are summarized in Table |, 
that shows the computed density p c for various values of 
TV and O along with the number of relevant configuration 
N c . An evident fact is that TV C grows very rapidly as a 
function of O and TV, making it impractical to perform 
the calculation for wide cylinders. We note that in order 
to obtain the same relative accuracy it is necessary to 
use O ex TV , e.g., in order to obtain a relative accuracy 
better than 10 -4 one should use at least O = TV — 1. This 
is the case for TV < 7, where the results are very accu- 
rate, but not so for TV > 8, where our available computer 
resources allo wed o nly lower order computations. As dis- 
cussed in Sec. [ID, we are able to improve the estimates 
in these cases by extrapolating O — > oo, taking advan- 
tage of the universal exponential decay of p c (N, O)fp(N) 
with the scaled variable O/N. Table | also compares the 
enumeration estimates with direct measurements from 
simulations, and finds them to agree within the simu- 
lation errors. Once accurate estimates are obtained for 
p(N) for TV < 10, they are fitted to a power-law approx- 
imati on w ith a correction to scaling term according to 
Eq. (2.31). The fit (with 8 = 1) gives an estimate of the 
fractal dimension D = 1.68 ± 0.01. 

Besides the range 2 < TV < 10, simulations are also 
performed on cylinders with larger TV's in the range 
10 < TV < 128. The relative errors of the measurements 
of p(N) are estimated around 10 -4 . The simulation data 
are also fitted to the same approximations. Once again, 
the three parameter approximation proves most appro- 
priate and the resulting fractal dimension this time is 
D = 1.671. The fact that the enumeration and simu- 
lation based estimates of the fractal dimension are very 
close is a good indication of their accuracy. 

The last statement should be taken with some caution 
in light of evidence that raises doubts concerning self- 
similarity in radial DLA Jl9|j20|j2l[| , or suggesting some 
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very slow crossovers 22 23j]. Indeed, radial DLA is some- 
what different than cylindrical DLA, as manifested by the 
difference between their fractal dimensions: D = 1.71 for 
radial DLA HIH and D = L66 for cylindrical DLA 
(this difference is still not fully understood). 

We also tried performing the exact calculations for 
TV = 11 and 12, but managed to go only up to O max = 4. 
This was insufficient for extrapolation with an accuracy 
that is comparable to the rest of the data points. With 
the aid of stronger computers we think that it would 
be possible and beneficial to compute a few more data 
points p(N), which would help obtaining more accurate 
estimates of the fractal dimension. Also, the techniques 
discussed here could be used to find the fractal dimen- 
sion of cylindrical DLA in 3D. However, since a much 
larger number of configurations can be expected, this 
task would also probably require the aid of a very strong 
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computer. 

There are a few differences between site-DLA and 
bond-DLA: The boundary conditions for the Laplace 
equation are a little bit different; In bond-DLA the poten- 
tial is set to zero on the aggregate itself, whereas in site- 
DLA the potential is set to zero on sites that are nearest 
neighbors of the aggregate. Also, the growth probabilities 
are computed somewhat differently; In bond-DLA con- 
tributions are summed over bonds that go out of a site 
where sticking may occur, whereas they are summed over 
bond that go into it in site-DLA. The normalization fac- 
tor, however, is equal to the width N is both cases. In the 
case of site sticking conditions there is an effective thick- 
ening of branches and thus a narrowing of fjords. Thus, 
there is a notable decrease in the probability of a random 
walker to penetrate deep into fjords. This also causes 
the number of configurations for a particular choice of N 
and O to be considerably less for site-DLA in comparison 
with bond-DLA. Therefore, accurate enumeration results 
can be obtained for larger iV's and O's in site-DLA. The 
extrapolation O — > oo performed in this paper was not 
done in Ref. pi, which deals with bond-DLA, because 
the technique was not developed at that time. When 
we apply the method to the bond-DLA case, we manage 
to improve the relative accuracy of the highest available 
approximations, p c (N, O max (N)) for N = 6, 7 by an or- 
der of magnitude: from about 1.2 x 10~ 3 to 2 x 10~ 4 
for N = 6, and from 5 x 10~ 2 to 1.6 x 10~ 3 . This ex- 
trapolation is based on the data points for N = 5. The 
relative accuracy of p c (N, O max (N)) for N < 5 is better 
than 10 -4 and hence, the extrapolation is not necessary. 
The estimate of the fractal dimension for site-DLA is, 
D = 1.68, to be compared with the bond-DLA enumer- 
ation result D = 1.64 |1J]. In contrast, the difference 
in the simulation results for the two cases is smaller: 
D = 1.67 for site-DLA and D = 1.66 for bond-DLA. 
Given the uncertainties, our results are consitent with 
universality with respect to the sticking conditions (lq,p[ . 
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FIG. 1. An example of the solution of the Laplace equation V 2 $(m, n) — 0, with boundary conditions $ = on the aggregate, 
and <9$/<9m = 1 on the upper distant boundary. Here, the width is N = 5 and there are periodic boundary conditions on the 
sides. The axes indicate the directions of the coordinates m and n. These boundary conditions are consistent with bond-sticking 
conditions. 
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FIG. 2. The growth probabilities for the aggregate shown in Fig. |l|. The growth probability in each perimeter site is 
proportional to the potential $ at that site and to the number of bonds JV& leading from the site into the aggregate (denoted 
by arrows), e.g., Nt = 3 for the site at ( — 1, 2) and Nt = 1 for the site (0, 2), right above it. 
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FIG. 3. The solution to the Laplace equation near the same aggregate as in Figs, hi andH, only with site-sticking conditions. 
The circles denote the perimeter sites where a random walker might stick. The boundary conditions are that $ = on these 
sites, unlike the case of bond-sticking conditions where $ = 0on the aggregate itself. The boundary conditions d^/dm — 1 at 
large m remain unchanged. 
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FIG. 4. Two different aggregates (represented by the dashed squares) with TV = 2, having exactly the same set of sites where 
a random walker may stick (shown in circles), and thus having the same boundary (the bold line) for the Laplace equation, 
where $ = 0. 
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FIG. 5. The sticking probabilities in each of the circled sites of Fig. pi They are computed by summing over all of the bonds 
that go into each site (denoted by arrows), unlike the case of bond-sticking conditions, where contributions are summed over 
bonds that go out of each site . 
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FIG. 6. The two possible configurations for N = 2. The circles denote the sites where a random walker might stick. Also 
shown are the possible transitions between them, denoted by arrows, and the relevant matrix elements Eij. The distribution 
of the potential $ over the lattice is demonstrated only for configuration i = 2, see text for explanation. The double arrows 
(•!)., => and 4=) show the bonds of the possible access paths, which a random walker can take into the circled sites. The bold 
double arrow shows the only bond going into site (0, 0). The other three bonds lead into site (1, 1). 
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FIG. 7. The possible configurations and the possible transitions between them for N — 3. 
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FIG. 8. A "potential diagram": the potentials <&{m,n) of configuration i — 2, expressed in terms of the variable x. 
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FIG. 9. The potential diagram for configuration i — 3. 
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FIG. 10. The only configurations for a cylinder of width N = 4: (a) Am = 2, (b) Am > 2. 
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FIG. 11. The four possible configurations for TV = 4 with Am = 1. These configurations are indexed between i — 2 and 
j = 5, and the flat configuration with Am = is indexed by i = 1. 
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FIG. 12. Data collapse of e f = \p c (N,0)/p(N) - 1] vs. O/N for all the data points with 4 < N < 10 and O > 2, on 
semi-log scale. The continuous line shows the linear approximation / ~ —12.30/N. 
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FIG. 13. The bold line separates between the upper inactive zone with m > 1 and the lower zone with m < 1. A random 
walker cannot stick in the inactive zone. 
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FIG. 14. The density profile as a function of height for a lattice with N = 10 averaged over some Ni ~ 2 x 10 7 iterations. 
The height of the lattice is 140 sites. 
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FIG. 15. A plot of p(N) vs. N on log- log scales. The plus sign s denote the simulation results, the dashed line denotes p a {N) 
the best fit to the three parameter approximation of Eq. (2.31), with 9 — 1. The circles denote the enumeration results. 
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FIG. 16. The relative residuals v — (p — pa)/ p (the plus signs) vs. N on semi-log scales. The upper and lower triangles show 
the estimated confidence intervals (errors) of the simulation data ,±cr. 
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TABLE I. The approximated densities p c and the number of configurations N c for various orders O and cylinder widths 
N. The approximated densities from enumeration are compared to simulation results. In addition, the extrapolated density 
p(N, O — ► oo ) is also presented. 
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